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Using density functional molecular dynamics free energy calculations, we show that the body-centered-cubic 
phase of superionic ice previously believed to be the only phase is in fact fhermodynamically unstable compared 
to a novel phase with oxygen positions in fee lattice sites. The novel phase has a lower proton mobility than 
the be phase and may exhibit a higher melting temperature. We predict a transition between the two phases 
at a pressure of 1 ± 0.5 Mbar, with potential consequences for the interiors of ice giants such as Uranus and 
Neptune. 
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Water is one of the most prevalent substances in the uni- 
verse and exists in a large number of phases over a vast range 
of temperature and pressure conditions. In addition to the liq- 
uid, gas, plasma and many solid phases [l-5 j, water also pos- 
sesses a superionic phase, in which the oxygen atoms occupy 
fixed lattice positions as in a solid, while hydrogen atoms mi- 
grate through the lattice as in a fluid Sill- The superionic 
phase is predicted to occupy a large section of the ice phase 
diagram for pressures in excess of 0.5 Mbar and temperatures 
of a few thousand Kelvin j6]-0. As this regime corresponds 
to conditions in the interiors of ice giants such as Uranus and 
Neptune, which are believed to consist largely of water, it is 
predicted that these planets consist largely of superionic ice 
Q, making an understanding of the physical and chemical 
properties of superionic ice vital for understanding the inte- 
rior structure and evolution of these planets. 

Although superionic ice has been extensively studied in ab 
initio theoretical studies (6]-|9], all works up to this point have 
assumed the superionic phase to maintain a body centered cu- 
bic (bec) structure for the oxygen sublattice, as seen in the 
solid ice VII and ice X phases [ 10 1. In this paper we predict 
instead, via density functional theory (DFT) free energy cal- 
culations, that the bec phase is thermodynamically unstable 
relative to a denser face centered cubic (fee) phase for pres- 
sures in excess of 1.0 ± 0.5 Mbar. The fee phase is found to 
have a lower hydrogen mobility than the bec phase. The pro- 
posed phase transition may intersect with the Neptunian and 
Uranuian isentropes, suggesting the possibility of a superionic 
to superionic phase transition in ice giants. 

The existence of superionic ice was initially predicted in 
via DFT-MD simulations by Cavazzoni et al [81, by heating 
of the ice X and ice VII phases of water at pressures in ex- 
cess of 0.5 Mbar. The bec oxygen sublattice of the ice X and 
VII was found to be maintained upon the transition to the su- 
perionic phase. Goldman f9j et al in 2005 studied bonding 
and diffusion in superionic water, again assuming the oxygen 
atoms to retain a bec sublattice. French et al. (6l EJ exten- 
sively studied the bec superionic phase and its transition to 
the fully fluid or plasma regime in which both hydrogens and 
oxygens become mobile. In repeated simulations, French et 



al. cooled water from the fully fluid regime and observed the 
re-formation of a superionic phase with a bec oxygen sublat- 
tice, however the geometric constraints of the unit cell used, 
with 54 H2O molecules in a cubic cell, mean that the forma- 
tion of alternative structures whose sublattices do not conform 
to these constraints is penalized. Experimentally, superionic 
ice has been observed in laser-heated diamond anvil cell ex- 
periments by Goncharov et al. ifTTl who demonstrated spec- 
troscopically a phase transition believed to correspond to the 
superionic phase at approximately 0.47 Mbar, however these 
experiments did not provide structural information. 

Hints of the instability of the bec oxygen sublattice over at 
least some portion of the superionic ice regime have been ob- 
served in several studies. French et al [6| noted the existence 
of a region of the phase diagram at low temperature and mod- 
erate pressure in which the bec oxygen sublattice was unstable 
within molecular dynamics (MD) - that is, at which the sys- 
tem readily evolves out of the sublattice in timescales acces- 
sible to our MD simulations. In the present authors' study of 
the solubility properties of superionic ice at giant planet core 
conditions [12] we briefly noted that the bec phase became un- 
stable at two sets of conditions under consideration (10 Mbar 
with temperatures of 2000 and 3000 K). Recent work on supe- 
rionic ammonia [ 13 ] has raised the possibility of the existence 
of phase changes within the superionic regime of NH3. These 
factors motivated a formal study of alternative oxygen sublat- 
tices in superionic water. 

The first stage of our study was to investigate the short- 
term stability of the bec lattice as a function of temperature 
and density, along with that of the fee and hexagonal close 
packed (hep) lattices. The ability of a system to retain a partic- 
ular structure over the picosecond timescales associated with 
a molecular dynamics simulation is a necessary but insuffi- 
cient condition for a structure to represent the thermodynamic 
ground state. We restricted our attention to these three sim- 
ple high-symmetry sublattices on the tentative assumption that 
thermal vibrations and the disordered migration of the hydro- 
gen atoms mean that the oxygen sublattice is likely to possess 
a configuration with high degree of symmetry. To investigate 
the short-term stability of a superionic lattice structure, we 
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first undertake a molecular dynamics simulation in which the 
oxygen atoms remain constrained to lattice positions while the 
hydrogen atoms move freely to equilibrate at the desired tem- 
perature. The constraint on the oxygen atoms is then released 
and the dynamics continued, with a newly-generated thermal 
velocity distribution. We simulated fee and bec ice struc- 
tures in constant cells at temperatures of 2000-5000 K and 
at nine densities from 3 gcm~ 3 (approximately 1 Mbar) to 1 1 
gem -3 (approximately 40 Mbar). We used the VASP package 
GTl and pseudopotentials of the projector-augmented wave 
type E2| with the exchange-correlation functional of Perdew, 
Burke and Ernzerhof ll23l . Supercells with 32, 54 and 48 
molecules were used for the fee, bec and hep structures re- 
spectively. Molecular dynamics simulations used a timestep 
of 0.2 fs and were integrated over between 2000 and 5000 
steps. A distortion in the bec oxygen sublattice was observed 
at 2000 K for densities of 6 and 7 gcm~ 3 (corresponding to 
pressures around 9 and 14 Mbar) and also at 6 gcm~ 3 at a 
temperature of 4000 K. In a later simulation we also found 
the bec structure to become distorted at a pressure of 40 Mbar 
and temperature of 10000K. The fee superionic structure re- 
mained stable in MD under all conditions studied. 

All attempts to simulate superionic ice with an hep oxy- 
gen sublattice rapidly resulted in the hep oxygen sublattice 
becoming distorted. Although the hep and fee lattices are very 
similar, differing only by the stacking of layers, and represent 
equivalent packings of spheres, we note that there is an im- 
portant difference in their distribution of interstitial sites. In 
the fee superionic structure we find that the hydrogen atoms 
largely concentrate around the tetrahedral interstitial sites, of 
which there are two per oxygen atom. Due to the different ar- 
rangement of oxygen atoms in the hep crystal, the equivalent 
tetrahedral interstitial sites are arranged in very closely spaced 
pairs, making the simultaneous occupation of all tetrahedral 
sites disfavored. 

Having established the short-term stability of the fee and 
bec phases across most of the superionic ice regime, we must 
now determine which phase has a lower Gibbs free energy 
G = U + PV — TS. We chose several representative points 
throughout the superionic regime, ranging from pressures of 1 
Mbar close to the onset of superionic behavior up to 40 Mbar 
corresponding to the approximate pressure of Jupiter's core; 
we chose points at which both bec and fee phases were stable, 
steering clear of the low-temperature regime around p = 6 
gem" 3 (« 10 Mbar). 

For the computation of Gibbs free energies we adopt a two- 
step coupling constant integration (CCI) method as recently 
applied by several authors lfl2l [T4l - [T6l and identical to that 
described in more detail for our earlier work on superionic 
water solubility [ 12 1. In this scheme the free energy of the sys- 
tem of interest is computed from the free energy of a simpler 
system whose free energy is known analytically via a thermo- 
dynamic integration in which the simpler system is gradually 
changed into the system of interest. For the analytic system 
to resemble the superionic phase we choose a system consist- 
ing of noninteracting harmonic oscillators at lattice sites for 
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0.180(11) 0.033(19) 
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TABLE I: Difference in free energy Gtcc — Gf cc between the fee 
and bec phases, and the PAV, AU and TAS components of the 
free energy difference. 



oxygen atoms, and a noninteracting ideal gas for the hydro- 
gen atoms. The difference in Helmholtz free energy between 
systems 1 and 2 governed by potential energy functions U\ 
and U2 is given by 

Fi= [ (U 1 -U 2 ) x dX + F 2 (1) 
Jo 

where the integral is taken over a set of systems governed 
by a linear combination of the forces from the two systems, 
U\ = (1 — A)?7i + XU 2 . For numerical efficiency, the inte- 
gration is taken in two steps: firstly from the full DFT system 
and to a system governed by a simple empirical potential cho- 
sen to closely resemble the dynamics of the DFT system, and 
then from the empirical system to the idealized noninteracting 
system. For our empirical potential we used a simple two- 
body spline-form potential generated by the force-matching 
methodology of Izvekov et al fl9l in combination with a har- 
monic potential anchoring oxygen atoms to their lattice sites. 
Appropriate volumes for the simulation cell were determined 
using the variable-cell methodology of Hernandez [20]. 

Results of the free energy calculations are shown in Table 
[I] We find the fee structure to have a lower Gibbs free energy 
under all conditions studied, with the exception of the 1 Mbar 
and 3000 K point where the energy difference between bec 
and fee is nearly zero. Figure [T] plots free energy differences 
between the two structures as a function of pressure at 3000 K 
and of temperature at 40 Mbar. A strong tendency towards 
greater stability of the denser fee structure as pressure is in- 
creased. We predict a transition from bee to fee stability at a 
pressure of 1.0 ± 0.5 Mbar. The error bar on the transition 
pressure means that the possibility that bec may have no sta- 
bility field at all is not excluded. Figure|2]shows the estimated 
location of the bec to fee transition in relation to the phase 
diagram. 

The free energy difference between two structures may be 
broken down into three components; an internal energy term 
AU, a volume term PAV and an entropic term —TAS. The 
resulting components are listed in Table 2. It is notable that the 
pressure-volume and internal-energy terms increasingly favor 
the fee structure as pressure is increased, while the entropic 
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FIG. 1: Gibbs free energy difference per molecule between the bcc 
and fee phases (solid line) shown as a function of pressure for a tem- 
perature of 3000 K, and as a function of temperature for a pressure 
of 40 Mbar. Also shown are the internal energy A(7 and pressure- 
volume PAV components of the free energy difference. 



term has the opposite sign. 

Given the different sizes of cells used for the bcc and fee 
phases and especially given the differing fc-point densities 
which result from using a 2 x 2 x 2 k-point grid for both cells, 
it is necessary to carefully check the convergence of the com- 
puted quantities with respect to k-point sampling and super- 
cell size. Doubling the size of the fee supercell in one direc- 
tion resulted in a change of less than 0.010 eV per molecule, 
and increasing the k-point sampling from 2x2x2 to 4x4x4 
changed the computed internal energy by less than 0.01 eV at 
p = 10. For the bcc supercell we found a change of less than 
0.006 eV/molecule in the internal energy when increasing the 
k-point mesh at p = 5 and a somewhat larger change ( 0.02 
eV/molecule) in the internal energy of the bcc structure at high 
densities of p = 10 (approximately 34 Mbar) due to the clo- 



FIG. 2: Phase diagram of water showing the superionic regime. The 
phase boundary from solid ice to superionic water indicated by a 
shading gradient indicating the degree of uncertainty due to the error 
bars. The triple point and melting line of bcc superionic ice are from 
Redmer et al (7), as are the indicated isentropes of Uranus and Nep- 
tune. The superionic regime is shaded to show the transition from 
bcc to fee stability at pressures of 1 ± 0.5 Mbar. The fec-superionic 
to fluid and solid-to-superionic lines were obtained by single-phase 
simulations in which the system was heated and cooled; we esti- 
mate a slightly higher melting temperature for the fee lattice than 
was found by Redmer et al for the bcc lattice. 



sure of the band gap in the bcc structure at high pressures. 

A key physical characteristic of the superionic phase is the 
hydrogen mobility, shown in Table II. We find hydrogen to 
diffuse more slowly in the fee than the bcc structure under 
all conditions. The greater packing density of the fee lattice 
allows fewer channels for hydrogens to diffuse from one site 
to another. Figure [3] shows isosurfaces of hydrogen density 
throughout molecular dynamics runs carried out at 40 Mbar 
and 5000 K. Apparent from these images is the fact that hy- 
drogen atoms in the fee structure are largely confined to tetra- 
hedral intersititial sites, while in the bcc structure may mi- 
grate more freely between two different types of interstitial 
site (tetrahedral and octohedral). The greater variety of sites 
available to the hydrogen atoms in the bcc structure may also 
explain the entropic preference for the bcc structure. The tran- 
sition from bcc to fee thus coincides with a drop in hydrogen 
mobility, with consequences for thermal and electrical con- 
ductivity properties. 

We have established the stability of the fee over the bcc su- 
perionic phase for pressures in excess of 1.0 ± 0.5 Mbar, im- 
plying a strong possibility of a superionic -to-superionic phase 
transition within this pressure regime. Such a transition could 
be experimentally observed by laser heated diamond anvil cell 
ifTTI or laser-driven shock experiments [17|. The rearrange- 
ment of oxygen atoms in a bcc to fee transition can best be de- 
tected with X-ray diffraction or X-ray Absorption Near Edge 
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FIG. 3: Isosurfaces of constant hydrogen density in molecular dy- 
namics simulations of superionic ice in the bcc (upper) and fee 
(lower) phases at 40 Mbar and 5000 K. Two different maxima repre- 
senting the tetrahedral and octohedral sites are seen in the bcc lattice; 
in the fee lattice by contrast the hydrogens are concentrated at the 
tetrahedral interstitial sites. 
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1.55 
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5.24 
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TABLE II: Diffusion constants for hydrogen in bcc and fee superi- 
onic ice, in A 2 fs _1 . 

Structure techniques. Figure [2] also implies that a superionic 
transition may occur along the isentropes of Uranus and Nep- 
tune, or that alternatively these planets may bypass the bcc 
superionic regime altogether. 

Previous theoretical studies of superionic water require 
some reconsideration in light of these results. The fee phase 
has a higher melting temperature than the bcc phase, which 
will change computed equations of state for this material and 



may lead to a larger superionic ice regime in giant planet in- 
teriors than had previously been considered ]6j 0. Interior 
models of Uranus and Neptune will require some revision, 
although the relatively small (< 0.5%) difference in density 
between the two phases may preclude a large effect. The 
consequences of a superionic to superionic phase transition 
should be considered in the context of whether such a transi- 
tion may help explain the non-axisymmetric non-dipolar mag- 
netic fields of these two planets ifTSl . The conclusions of our 
previous work on the solubility of water ice in metallic hydro- 
gen within gas giant planet interiors [12| do not change sig- 
nificantly, due to the small magnitude of the free energy dif- 
ference between superionic phases (0.1 eV) compared to the 
large free energy associated with solubility at Jupiter and Sat- 
urn core-mantle boundary conditions. Further work, includ- 
ing understanding the potential implications of phase changes 
in the superionic regime for the convective and heat trans- 
port properties of Uranus and Neptune, as well as experi- 
mental work aimed at detecting this phase change in practice, 
may provide further insight into the interiors of these poorly- 
understood ice giants. 
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